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ABSTRACT 

Digital enhancement of cometary images is a necessary tool in studying cometary morphology. 
Many image processing algorithms, some developed specifically for comets, have been used to 
enhance the subtle, low contrast coma and tail features. We compare some of the most commonly 
used algorithms on two different images to evaluate their strong and weak points, and conclude 
that there currently exists no single “ideal” algorithm, although the radial gradient spatial filter 
gives the best overall result. This comparison should aid users in selecting the best algorithm to 
enhance particular features of interest. 


INTRODUCTION 

The observed morphology of cometary comae and tails is determined by ejection circumstances 
and the interaction of the ejected material with the cometary environment. Temporal data on the 
anisotropic emission of dust and gas over time can provide useful information on such things as 
location of active areas on the nucleus, the nucleus spin vector, and gas/dust interaction. Time 
sequences may also be particularly useful in following the formation and evolution of fast moving 
ion structures near the nucleus. However, discrete coma features are usually diffuse, of low 
amplitude, and they are superimposed on a steep intensity gradient radial to the nucleus. The 
usual 8-bit display does not usually show much detail because for large intensity ranges the 
contrast of coma structure is too small, or if the contrast is increased, only a small intensity range 
is visible. 

To improve the visibility of these features, a variety of digital enhancement algorithms have been 
employed with varying degrees of success. They usually produced some spatial filtering, and 
were chosen to optimize visibility of certain detail. Since information in the image is altered, it is 
important to understand the effects that parameter selection and processing artifacts can have on 
subsequent interpretation. Our definition of the ideal algorithm is that it must enhance low 
contrast features while not introducing misleading artifacts. We require that features seen in the 
enhanced image also be seen in the unenhanced version with appropriate contrast and intensity 
windowing. We have not found the ideal algorithm; all of those that we have tried or developed 
have strong and weak elements depending upon the type of features of interest. To help assess the 
suitability of various algorithms, we have processed two images in a variety of ways so they can 
be directly compared. We attempt to identify the strong and weak points of each in the context of 
optimizing visibility while maintaining positional integrity of features at the expense of 
photometric information. 

Figures 1 and 2 show the effects of various image processing on two sample images containing a 
variety of coma and tail features. An attempt has been made to generate hard copy that shows as 
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much detail as possible for each image. We do not assert that these results will be the most 
appropriate for other data. Experience with many data sets taken with a variety of telescopes and 
plate scales indicate that the sample frequency (pixel size) to resolution ratio can be important in 
selecting the algorithm or processing parameters. In general, there is more parameter leeway with 
over-sampled data (pixel size/”resolution”< 1) than for under-sampled data because of reduced 
centroid Cor registration) uncertainty and pixelization artifacts are smaller relative to the features 
of interest. The two sample images and their processed versions are shown in Figures 1 A-H and 2 
A-H. The figure letters correspond to the subheading letters below. 

A. THE UNPROCESSED IMAGES 

The sample image in Figure 1 is a composite which includes a mostly H20 + CCD image of P / 
Brorsen-Metcalf taken with a narrow-band filter transmitting the H20 + (0,8,0) band near 619 nm, 
and some of the C 2 (0,2) band. This comet was extremely dust-poor, so the ion streamers were 
relatively high contrast, and the C 2 is easily seen to be restricted near the nucleus. However, since 
this is not usually the case for a “typical” dusty comet, we have added a 1/r intensity gradient to 
the central condensation that might be expected from isotropic emission of dust from a rapidly 
rotating nucleus. The resulting test image has a “dust”/ion intensity ratio of 50:1. A second 
combination using an image of P/Brorsen-Metcalf taken 15 minut es later was used for the 
temporal derivative. The scale of the image is 0.7 arcsec per pixel, or 400 km per pixel at the 
comet. 

The sample image in Figure 2 is a broad red-band CCD image of P/Halley that contains several 
dust jets, an antisunward dust jet, and ion streamers. A second image taken 45 minutes later was 
used for the temporal derivative. Both are three co-added 120 second exposures, and are also 0.7 
arcsec per pixel, or 640 km per pixel at the comet. 

B. LOG 10 INTENSITIES 

The simplest operation is to alter the intensity scale by some non-linear contrast stretch (such as 
the base 10 logarithm) to suppresses the steep intensity peak near the photocenter. In this case, the 
resulting pixel value (I) becomes; 

I (x,y) = 1 °g I o(x,y) 

where L, = original pixel count value. This helps bring the brightness of the faint outer regions of 
the coma closer to those in the inner coma. No artifacts are introduced, but there is only a mode st 
improvement in visibility of low contrast features. In gen era l, such non-lin ear str etches do not 
allow appreciable increase in local feature contrast relative to die background. The subjective 
result is to convert linear-response electronic images to a photographic-like response. The Strong 
point is that the image is easily understandable and interpretation is straight forward. Because the 
operation involves no special parameter selection (or “interpretation”) , t his is the form Near 
Nucleus Studies Network images appear in the International Halley Watch printed archive. 

SPATIAL DERIVATIVES 


C. Linear Shift-difference 

Spatial derivative, or “shift-difference” algorithms enhance intensity discontinuities, but only in 
the direction of the shift (Klinglesmith, 1982); 

I(x, y) = ^(x, y) " ^o(x+n, y+m) 

where, n, m are the amount of the shift in x and y. The simple linear shift-difference is very 
effective for bringing out ion tail structure when the shift is perpendicular to the tail axis. The 
magnitude of shift is a compromise between showing the smallest detail and enhancing the noise. 
The most serious problems with this method is that the results are directionally dependent, and are 
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not easily interpretable. The resulting features, showing the rate of change of intensity, enhance 
the edges of jets and shells. In both of these examples, the shift is perpendicular to the tail axis to 
best show the ion features. The residual radial gradient is directly related to the magnitude of 
shift. However, a smaller shift would leave only the highest spatial frequencies (mostly system 
noise). Coma detail is particularly difficult to interpret. 

D. Rotational Shift-difference 

Radial and rotational shift-difference algorithms were developed in response to the fact that 
particles generally move radially away from the nucleus, and that nucleus rotation produces spiral 
features (Larson and Sekanina, 1984). In polar coordinates, the combination radial and rotational 
shift gives; 

fy.e) = W.O) ' X ((r+Ar), (0+A0)) 

where Ar = the radial shift, and A0 = the rotational shift. The rotational derivative used in this 
comparison emphasizes rotational intensity discontinuities, such as radial jets and ion streamers. 
A radial derivative emphasizes discontinuities concentric to the nucleus, such as hoods and 
evolved jets. In one variation, such shifts are done in combination with rotational shifts of differ- 
ent directions but different amounts to reduce the directional dependency of the enhancement. 
The reproduced images show features in a 10° rotational derivative, however, ion streamer over- 
lap may cause confusing aliasing. The angular shift produces a radial dependency on the linear 
scale of edges that become enhanced. For some types of gas jets that become more diffuse farther 
from the nucleus, rotational shift-differencing is useful, but for ion features, it is usually confus- 
ing. With the rotational derivative, it is possible to see features down to the central condensation. 

Since the synthetic test image contains no radial discontinuities, nothing appears in a radial shift- 
difference, but the Comet Halley image shows the outer envelope well. By combining the radial 
plus rotational shift difference the directional dependency of feature visibility is reduced but not 
eliminated. Interpretation requires great care since it is the edges of features that are enhanced. 
Finding the best combination of rotational and radial shifts can be difficult and enhances only a 
limited spatial frequency range. Radial and combination derivatives are not shown in the figures. 

E. TEMPORAL DERIVATIVE 

Differencing carefully registered images taken at different times yield features that have moved. 
Such temporal derivative images are projected velocity maps that (among other things) make it 
easy to distinguish rapidly varying ion features from the slower moving dust structures (e.g. 
Larson and Minton, 1971; Larson, 1986). 

X (x, y) = Jofxtj, ytj) • I o(xt 2 . y^) 

where tj and t 2 are different times. Successful short-term difference images bring out ion features 
which are normally a minor component of a broad-band image. This method places great 
demands upon sets of images with nearly identical quality and very precise registration. 
Registration on the centroid of a several pixel area around the central condensation is usually 
necessary to reduce the effects of noise. Variations in seeing and guiding can complicate the 
result. Interpretation must be made with care, since the result is an image of moving feature edges. 

The synthetic test image processed this way shows the ion features best of all of the algorithms 
used, as all of the radial intensity gradient is eliminated. The rapid ion motions require images 
taken only a few minutes apart to show detail at this resolution. The Comet Halley image also 
shows the ion streamer motions, but also a lot of garbage resulting from imperfect flat fields, 
scattered twilight, and possible differential extinction at the large airmass. These problems cause 
the coma to have dark and light components that have nothing to do with the comet. The best 
temporal derivatives resulting from consistent, regularly spaced sequences of images are 
relatively easy to interpret 
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F. AZIMUTHAL RENORMALIZATION 

Azimuthal function algorithms reduce the radial gradient by either subtracting the average value 
in the annulus of constant distance from the photocenter, or by subtracting a best-fit, low- order 
function to the annulus (A’Heam, et al., 1986). 

I(r, 6) = Io(r, 0) - ((Z I(r))/n)) 

where n = number of pixels falling within die annulus of radius r. and are summed overall theta. 
Averaging and function fitting is more easily done after a polar to rectangular coordinate 
transformation (with the photocenter at the origin). This method is very efficient in eliminating 
the radial gradient, does not have any directional dependencies, and interpretation is straight- 
forward. The photocenter must be determined very carefully, especially for undersampled data, or 
spurious features close to the nucleus may be produced. Options in the A’Heam and Klavetter 
(Univ. Maryland) program includes normalization to mean or median values at constant radius, or 
normalizing to a functional profile. Our own implementation of this algorithm normalizes to the 
mean of the annulus without a polar to rectangular coordinate conversion step (fig. 2F). 

The synthetic image responds very well by eliminating all of the radial gradient and leaving the 
ion features intact. The noise becomes rather apparent, but there are options for averaging over 
various ranges of r and theta that the authors did not have time to try. We ar e not confident that the 
optimum parameters for this algorithm were used, however, the result is the most easily 
interpretable version for this image. The Comet Halley image, primarily because of the large 
amplitude of the sunward asymmetry, does not respond as well to an azimuthal normalization to a 
mean value. Similarly, normalization to an azimuthally symmetric function produces over- 
compensation in the tail direction, and under-compensation in the sunward direction. Although 
this algorithm does not work very well for primarily continuum images it works very well on 
images of the more symmetrical gas comets. 

Another approach is to subtract a synthetic image based on a generalize d model of particle . .. 
outflow. TTiis assumes some a priori knowledge of the ejection function. If there is enough data, 
subtraction of a mask produced by the median of many images over time may also work (Jewitt, 
1991). In a sense, this is similar to a temporal derivative. 

G. SPATIAL FILTERING 

Traditional spatial filtering algorithms reduce the radial gradient by eliminating the low spatial 
frequency domain in the image. A “high-pass” gaussian deconvolution retains features in the 
image smaller than the gaussian, and by eliminating the broad radial gradient, the contrast of the 
smaller features can be increased. There is no directional dependency, but “ringing” artifacts can 
be seen around bright stars and the central condensation of the comet. Selecting the optimum size 
gaussian usually depends upon the characteristic size of the features of interest. Although large, 
complete filter kernels are more rigorous, our spatial filtering algorithm uses an abbreviated^ 
kernel convolution to reduce the computation time. A 3 x 3 pixel kernel is used, but the kernels 
size is adjusted by varying the spacing between kernel pixels. This allows filtering over a wide 
range of spatial domains, while operating on a small number of pixels. The kernel shape is 
dictated by weights given the “comer” and “medial” pixels in the kernel. 

I(r, 0) = Io(r,e) - SI o(r\9’) w l ' SI o(r”,0”) w 2 

where; w^ = “comer” weights; W 2 = “medial” weights, nominally (wj + W 2 ) = 1.0; 

Ax = Ay = kernel pixel separation, and r = (x 2 + y 2 ) 1 ^, r’ = ((x ± Ax) 2 + (y + Ay) 2 ) 1 ^ 2 ; 
r” = (x 2 + (y + Ay) 2 ) 1/2 ; ((x ± Ax) 2 + y 2 ) 1 ^ 2 , 8 = arctan (y/x); 0’ = arctan ((y + Ay)/(x ± Ax)), 

8” = arctan (y/(x ± Ax)), ((y + Ay)/x) 

Processing the synthetic image brings out the ion tail Fairly well up to the point that the spatial 
component of the radial gradient becomes smaller than the effective kernel FWHM near the" 
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central condensation. Use of a smaller kernel to reduce this residual gradient emphasizes the noise 
and is better suited to high S/N ratio images. The P/Halley image shows nearly all of the jet, 
envelope and ion features quite well. The kernel size was chosen to optimize over-all visibility of 
features, but there is a residual small gradient close to the central condensation. Interpretation of 
this image is not difficult. 

H. RADIAL GRADIENT SPATIAL FILTERING 

Since cometary features typically become larger farther from the nucleus, spatially selective 
spatial filtering might be desirable. We have recently developed a variable kernel deconvolution 
routine that passes increasingly higher frequencies closer to the photocenter. This enhances a 
larger spatial range of coma features which often exist in an image. The previous algorithm is 
used, but the spacing (Ax, Ay) is variable from the photocenter, 

Ax = Ay = (a + (r/b)), a = kernel pixel separation, b = radial scale factor. 

The processed synthetic image shows the ion tail well, but does not show the folding rays as well 
as the temporal derivative. Note the strong pixelization pattern radiating from the nucleus. The 
processed P/Halley image shows detail of the central condensation as well as the ion streamers. 
The parameters wj, W 2 , a and b must be chosen empirically to ensure that features in the various 
spatial regimes are visible. Field stars may present problems with “ringing” (dot patterns in this 
abbreviated kernel). 


CONCLUSIONS 

We find that for coma jets, spatial filtering provides the best overall enhancement with minimal 
artifacts, and the results are readily interpretable. The radially variable convolution kernel 
provides a range of spatial filtering appropriate for different regions of the comet (high-pass at the 
center, lower near the edges). For ion structure, we find that temporal derivatives are the most 
useful in not only bringing out low contrast detail, but providing a “velocity map” of ion motions. 
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Fig. 1. A. Test image with ion features and a symmetric 1/r intensity gradient (see text), B. 
displayed as the base 10 logarithm of the counts, C.linear derivative (perpendicular to the tail 
axis), D. rotational derivative, E. temporal derivative, F. azimuthal renormalization, G. constant 
spatial filter, and H. radially varying spatial filter. 
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Fig. 2 A. Image of P/Halley (see text), B. displayed as the base 10 logarithm of the counts, 
C.linear derivative (perpendicular to the tail axis), D. rotational derivative, E. temporal derivative, 
F. azimuthal renormalization, G. constant spatial filter, and H. radially varying spatial filter. 



